################################################################################
# Replication File for
# "Party Elites' Preferences in Candidates: Evidence from a Conjoint Experiment" 
# Author: Jochen Rehmert
# jochen.rehmert@uzh.ch
################################################################################


################################################################################
#                              preparation                                
################################################################################


setwd('')

pkgs <- c("cregg","grid","ggplot2","readstata13","multiwayvcov","lmtest",
          "stargazer","MASS","cowplot")

## Install uninstalled packages 
# WARNING: This will install packages on your machine
lapply(pkgs[!(pkgs %in% installed.packages())], install.packages)

## Load all packages to library and adjust options
lapply(pkgs, library, character.only = TRUE)


### Data preparation
dat <- read.dta13("conjoint_data.dta")

dat$age[!is.na(dat$age) & dat$age ==1] <- "30 Years Old"
dat$age[!is.na(dat$age) & dat$age ==2] <- "45 Years Old"
dat$age[!is.na(dat$age) & dat$age ==3] <- "65 Years Old"

for(cc in 4:13){
  dat[,cc] <- as.character(dat[,cc])
}

# translate attribute levels into English
dat$gender[dat$gender == "Mann"] <- "Male"
dat$gender[dat$gender == "Frau"] <- "Female"
dat$education[dat$education == "Hauptschule"] <- "Secondary School"
dat$education[dat$education == "Abitur"] <- "High School Diploma"
dat$occupation[dat$occupation == "Technischer Beruf"] <- "Technical Profession"
dat$occupation[dat$occupation == "Rechtswesen"] <- "Legal Profession"
dat$occupation[dat$occupation == "Lehrer"] <- "Teacher"
dat$occupation[dat$occupation == "Landwirt"] <- "Farmer"
dat$district[dat$district == "Nein, keine Wahlkreiskandidatur"] <- "No District Nomination"
dat$district[dat$district == "Ja, kandidiert in einem Wahlkreis"] <- "District Nomination"
dat$experience[dat$experience == "MdL"] <- "Member of State Legislature"
dat$experience[dat$experience == "Keine"] <- "None"
dat$experience[dat$experience == "Kreistag"] <- "Municipal Council"
dat$party_length[dat$party_length == "7 Monate"] <- "7 Months"
dat$party_length[dat$party_length == "6 Jahre"] <- "6 Years"
dat$party_length[dat$party_length == "2 Jahre"] <- "2 Years"
dat$ideology[dat$ideology == "verortet sich links in der Partei"] <- "Leftist Member of Party"
dat$ideology[dat$ideology == "verortet sich rechts in der Partei"] <- "Rightist Member of Party"
dat$ideology[dat$ideology == "verortet sich in der Mitte der Partei"] <- "Centrist Member of Party"
dat$party_office[dat$party_office == "Kein Parteiamt"] <- "No Office"
dat$party_office[dat$party_office == "Kreisverbandsvorsitz"] <- "Party Chair, Municipality"
dat$party_office[dat$party_office == "Ortsverbandsvorsitz"] <- "Party Chair, Village"
dat$party_office[dat$party_office == "Landesvorstand"] <- "Member of Regional Board"
dat$activity[dat$activity == "selten"] <- "Rarely Active"
dat$activity[dat$activity == "h�ufig"] <- "Often Active"

for(cc in 4:13){
  dat[,cc] <- as.factor(dat[,cc])
}


dat <- dat[!is.na(dat$party_id),]

dat$party_id <- as.character(dat$party_id)

# reorder factor levels
dat$age <- factor(dat$age, levels = c("65 Years Old", "45 Years Old", "30 Years Old"))
dat$party_length <- factor(dat$party_length, levels = c("6 Years", "2 Years", "7 Months"))
dat$ideology <- factor(dat$ideology, levels = c("Rightist Member of Party", "Centrist Member of Party", "Leftist Member of Party"))
dat$party_office <- factor(dat$party_office, levels = c("Member of Regional Board","Party Chair, Municipality", "Party Chair, Village", "No Office"))



### weights
# delegate numbers: CDU (64/499), SPD (71/456), FDP (30/600), Die LINKE (73/385), B90/Gr (100/330)
# total: 2270
# CDU: 7.796875 (inverse of 64/499)
# SPD: 6.422535
# FDP: 20
# LINKE: 5.273973
# B90/Gr: 3.3

dat$freqw <- NA
dat$freqw[dat$party == "CDU"] <- 1/(64/499)
dat$freqw[dat$party == "SPD"] <- 1/(71/456)
dat$freqw[dat$party == "FDP"] <- 1/(30/600)
dat$freqw[dat$party == "Die LINKE"] <- 1/(73/385)
dat$freqw[dat$party == "B�ndnis 90/Die Gr�nen"] <- 1/(100/330)



###############################################################################
#                                Figure 1                                     #
###############################################################################
all <- as.formula(paste("choice ~ age + gender + education * occupation + 
                        district + experience * party_length + ideology +  party_office + activity"))

### pooled model
mm.all <- mm(dat, all, id = ~ResponseID,  alpha = 0.05)

ggData <- data.frame(mm.all)

ggData$feature <- as.character(ggData$feature)

ggData$level <- factor(as.character(ggData$level), 
                       levels = rev(c("30 Years Old","45 Years Old","65 Years Old",
                                      "Female","Male", "Secondary School","High School Diploma",
                                       "Technical Profession","Legal Profession","Teacher","Farmer",
                                       "District Nomination","No District Nomination",
                                       "Member of State Legislature","Municipal Council","None",
                                       "7 Months","2 Years", "6 Years",
                                       "Leftist Member of Party","Rightist Member of Party",
                                        "Centrist Member of Party", "No Office","Party Chair, Village",
                                        "Party Chair, Municipality","Member of Regional Board",
                                         "Rarely Active","Often Active")))


ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "Pooled Sample"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData$feature2 <- NA
ggData$feature2[ggData$feature == "district"] = "Nomination"
ggData$feature2[ggData$feature == "experience"] = "Experience"
ggData$feature2[ggData$feature == "party_length"] = "Membership"
ggData$feature2[ggData$feature == "party_office"] = "Party Office"
ggData$feature2[ggData$feature == "activity"] = "Activity"

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

# plotting
p = ggplot(ggData[ggData$feature %in% c("district","experience","party_length","party_office","activity"),])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(feature2 ~ stat_label, scales = "free") + xlab("") + ylab("")#+ ylab("Marginal Means") 
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p

###############################################################################
#                                Figure 2                                     #
###############################################################################
dat$veteran <- ifelse(dat$time_delegate >= 2, 1, 0)

mm.first <- mm(dat[!is.na(dat$veteran) & dat$veteran == 0,], all, id = ~ ResponseID, alpha = 0.1)

mm.veteran <- mm(dat[!is.na(dat$veteran) & dat$veteran == 1,], all, id = ~ResponseID, alpha = 0.1)


ggData <- data.frame(rbind(cbind(fold = "First-time,\n n = 133", mm.first),
                           cbind(fold = "Veteran,\n n = 116", mm.veteran)))

ggData$level <- factor(as.character(ggData$level), 
                       levels = rev(c("30 Years Old","45 Years Old","65 Years Old",
                                    "Female","Male", "Secondary School","High School Diploma",
                                     "Technical Profession","Legal Profession","Teacher","Farmer",
                                     "District Nomination","No District Nomination",
                                      "Member of State Legislature","Municipal Council","None",
                                      "7 Months","2 Years", "6 Years",
                                      "Leftist Member of Party","Rightist Member of Party",
                                      "Centrist Member of Party", "No Office","Party Chair, Village",
                                      "Party Chair, Municipality","Member of Regional Board",
                                      "Rarely Active","Often Active")))
ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData$feature2 <- NA
ggData$feature2[ggData$feature == "district"] = "Nomination"
ggData$feature2[ggData$feature == "experience"] = "Experience"
ggData$feature2[ggData$feature == "party_length"] = "Membership"
ggData$feature2[ggData$feature == "party_office"] = "Party Office"
ggData$feature2[ggData$feature == "activity"] = "Activity"


ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

p = ggplot(ggData[ ggData$feature %in% c("district","experience",
                                  "party_length","party_office","activity"),])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(feature2 ~ fold, scales = "free_y") + xlab("") + ylab("")#+ ylab("Marginal Means") 
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p



###############################################################################
#                                Figure 4                                     #
###############################################################################
dat$num_valence <- 0

dat$num_valence <- ifelse(dat$district == "District Nomination", dat$num_valence + 1, dat$num_valence)
dat$num_valence <- ifelse(dat$experience == "Member of State Legislature", dat$num_valence + 1, dat$num_valence)
dat$num_valence <- ifelse(dat$party_length == "6 Years", dat$num_valence + 1, dat$num_valence)
dat$num_valence <- ifelse(dat$party_office == "Member of Regional Board", dat$num_valence + 1, dat$num_valence)
dat$num_valence <- ifelse(dat$activity == "Often Active", dat$num_valence + 1, dat$num_valence)
  
dat$del_gender[dat$del_gender == "Man"] <- "Male"
dat$del_gender[dat$del_gender == "Woman"] <- "Female"
dat$del_ideology <- NA
dat$del_ideology <- dat$lr_own - dat$lr_party

dat$num_homophily <- 0
dat$num_homophily <- ifelse(dat$gender == dat$del_gender, dat$num_homophily + 1,dat$num_homophily)
dat$num_homophily <- ifelse(dat$education == "High School Diploma" & dat$del_education == "studied", dat$num_homophily + 1,dat$num_homophily)
dat$num_homophily <- ifelse(dat$age == "30 Years Old" & dat$del_age <= 37.5  , dat$num_homophily + 1,dat$num_homophily)
dat$num_homophily <- ifelse(dat$age == "45 Years Old" & dat$del_age >= 37.5 & dat$del_age < 55 , dat$num_homophily + 1,dat$num_homophily)
dat$num_homophily <- ifelse(dat$age == "65 Years Old" & dat$del_age >= 55 , dat$num_homophily + 1,dat$num_homophily)
dat$num_homophily <- ifelse(dat$ideology == "Rightist Member of Party" & dat$del_ideology > 0 , dat$num_homophily + 1,dat$num_homophily)
dat$num_homophily <- ifelse(dat$ideology == "Leftist Member of Party" & dat$del_ideology < 0 , dat$num_homophily + 1,dat$num_homophily)


summary(out.hom <- lm(choice ~  veteran*num_homophily + as.factor(ResponseID) + as.factor(TaskID), data = dat))

summary(out.val <- lm(choice ~  veteran*num_valence + as.factor(ResponseID) + as.factor(TaskID), data = dat[row.names(out.hom[["model"]])  , ]))

out.val$vcov <- cluster.vcov(out.val, ~ ResponseID)
out.val.se <- coeftest(out.val, out.val$vcov)[,2]

out.hom$vcov <- cluster.vcov(out.hom, ~ ResponseID)
out.hom.se <- coeftest(out.hom, out.hom$vcov)[,2]

### Table 3
stargazer(list(out.val, out.hom), 
          se = list(out.val.se, out.hom.se), 
          omit = c("TaskID","ResponseID"))

### plotting
seeds <- 9
set.seed(seeds)
b.sim <- mvrnorm(1000, coef(out.val)[c("(Intercept)","veteran","num_valence","veteran:num_valence")], out.val$vcov[c("(Intercept)","veteran","num_valence","veteran:num_valence"),c("(Intercept)","veteran","num_valence","veteran:num_valence")] )

df.0 <- data.frame(1, 
                   veteran = 0,
                   num_valence = 0:5,
                   veteran_valence = 0 * 0:5)

df.1 <- data.frame(1, 
                   veteran = 1,
                   num_valence = 0:5,
                   veteran_valence = 1 * 0:5)
preds0 <- as.matrix(df.0) %*% t(b.sim)
preds1 <- as.matrix(df.1) %*% t(b.sim)

ame <- apply(preds1 - preds0, 1, mean)
cis95 <- apply(preds1 - preds0, 1, quantile, c(0.025, 0.975))
cis90 <- apply(preds1 - preds0, 1, quantile, c(0.05, 0.95))
 
par(mfrow = c(1,2), mai = c(1,.8,.2,.2),
    oma=c(1.5,1.5,1.5,.5))
plot(ame ~ c(0:5), frame = FALSE, xaxt = "n", yaxt = "n", pch = 19, ylim = c(-0.25,0.3),
     ylab = "Marginal Effect", xlab = "Number of Valence Cues", cex.lab = 1.5)
segments(c(0:5) , cis95[1,], c(0:5), cis95[2,], lty = 2)
segments(c(0:5) , cis90[1,], c(0:5), cis90[2,])
lines(rep(0,7) ~ c(-.5:5.5), col ="red", lty = 2)
set.seed(1)
rug(jitter(out.val[["model"]]$num_valence))
text(x = c(0:5), par("usr")[3], labels = c(0:5),  pos = 1, xpd = TRUE, cex = 1.5)
text(y = c(-0.2,0,0.2), par("usr")[3], labels = c(-0.2,0,.2),  pos = 2, xpd = TRUE, cex = 1.5)


set.seed(seeds)
b.sim <- mvrnorm(1000, coef(out.hom)[c("(Intercept)","veteran","num_homophily","veteran:num_homophily")], out.hom$vcov[c("(Intercept)","veteran","num_homophily","veteran:num_homophily"),c("(Intercept)","veteran","num_homophily","veteran:num_homophily")] )

df.0 <- data.frame(1, 
                   veteran = 0,
                   num_homophily = 0:4,
                   veteran_homophily = 0 * 0:4)

df.1 <- data.frame(1, 
                   veteran = 1,
                   num_homophily = 0:4,
                   veteran_homophily = 1 * 0:4)
preds0 <- as.matrix(df.0) %*% t(b.sim)
preds1 <- as.matrix(df.1) %*% t(b.sim)

ame <- apply(preds1 - preds0, 1, mean)
cis95 <- apply(preds1 - preds0, 1, quantile, c(0.025, 0.975))
cis90 <- apply(preds1 - preds0, 1, quantile, c(0.05, 0.95))

plot(ame ~ c(0:4), frame = FALSE, xaxt = "n", yaxt = "n", pch = 19, ylim = c(-0.25,0.3),
     ylab = "", xlab = "Number of Homophilous Traits", cex.lab = 1.5)
segments(c(0:4) , cis95[1,], c(0:4), cis95[2,], lty = 2)
segments(c(0:4) , cis90[1,], c(0:4), cis90[2,])
lines(rep(0,7) ~ c(-.5:5.5), col ="red", lty = 2)
set.seed(1)
rug(jitter(out.val[["model"]]$num_valence))
text(x = c(0:4), par("usr")[3], labels = c(0:4),  pos = 1, xpd = TRUE, cex = 1.5)



###############################################################################
#                               Figure 3                                      #
###############################################################################

### Gender 
mm.woman <- mm(dat[dat$del_gender == "Female"   ,], all, id = ~ResponseID)
mm.man <- mm(dat[dat$del_gender == "Male"  ,], all, id = ~ResponseID)


ggData <- data.frame(rbind(cbind(fold = "   Female,    \n n = 119", mm.woman),
                           cbind(fold = "     Male,    \n n = 131", mm.man)))

ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData$feature2 <- NA
ggData$feature2[ggData$feature == "gender"] = "Gender"

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0


p = ggplot(ggData[ ggData$feature %in% c("gender"),])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(feature2 ~ fold, scales = "free") + xlab("") + ylab("")
p.gender = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p.gender


### Age

mm.old <- mm(dat[!is.na(dat$del_age) & dat$del_age >= mean(dat$del_age  , na.rm = TRUE),], all, id = ~ResponseID)

mm.young <- mm(dat[!is.na(dat$del_age) & dat$del_age <= mean(dat$del_age  , na.rm = TRUE),], all, id = ~ResponseID)

ggData <- data.frame(rbind(cbind(fold = "       Old,      \n n = 104", mm.old),
                           cbind(fold = "    Young,    \n n = 116", mm.young)))

ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData$feature2 <- NA
ggData$feature2[ggData$feature == "age"] = "Age"

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0


p = ggplot(ggData[ ggData$feature %in% c("age"),])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(feature2 ~ fold, scales = "free") + xlab("") + ylab("") 
p.age = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p.age

### Education 
mm.high <- mm(dat[dat$del_education == "studied",], all, id = ~ResponseID)

mm.low <- mm(dat[dat$del_education == "not studied",], all, id = ~ResponseID)

ggData <- data.frame(rbind(cbind(fold = "High-Educated,\n n = 124", mm.high),
                           cbind(fold = "Low-Educated, \n n = 90", mm.low)))

ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData$feature2 <- NA
ggData$feature2[ggData$feature == "education"] = "Education"

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

p = ggplot(ggData[ ggData$feature %in% c("education"),])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(feature2 ~ fold, scales = "free") + xlab("") + ylab("")#+ ylab("Marginal Means") 
p.education = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p.education


### Ideology 

mm.right <- mm(dat[!is.na(dat$del_ideology) & dat$del_ideology > 0,], all, id = ~ResponseID)
mm.left <- mm(dat[!is.na(dat$del_ideology) & dat$del_ideology < 0,], all, id = ~ResponseID)

ggData <- data.frame(rbind(cbind(fold = "Right-leaning,\n n = 77", mm.right),
                           cbind(fold = "Left-leaning, \n n = 105", mm.left)))


ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData$feature2 <- NA
ggData$feature2[ggData$feature == "ideology"] = "Ideology"

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

p = ggplot(ggData[ ggData$feature %in% c("ideology"),])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(feature2 ~ fold, scales = "free") + xlab("") + ylab("")#+ ylab("Marginal Means") 
p.ideology = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p.ideology

# ----- Combine Graphs ----- # 

plot_grid(p.gender, p.age, p.education, p.ideology, ncol = 2, align= "v")



###############################################################################
#                                 Table 2                                     #
###############################################################################
IDs <- unique(dat$ResponseID)
dat2 <- dat[which(dat$ResponseID %in% IDs & dat$TaskID == 1 & dat$ConceptID == 1),]

summary(vet.0 <- glm(veteran ~ party + del_gender  + del_age , data = dat2[dat2$del_gender != "",],
            family = binomial("logit")))

summary(vet.1 <- glm(veteran ~ party + del_gender  + del_age +del_education, data = dat2[dat2$del_education != "" & dat2$del_gender != "",],
                     family = binomial("logit")))


summary(vet.3 <- glm(veteran ~ party + del_gender  + del_age + del_exp_dichotom, data = dat2[dat2$del_gender != "",],
                     family = binomial("logit")))

dat2$del_pty_any <- ifelse(dat2$del_ptyoffice != "" & dat2$del_ptyoffice %in% c("federal","regional"),1,0 )
dat2$del_pty_none <- ifelse(dat2$del_ptyoffice == "none",1,0 )

summary(vet.4 <- glm(veteran ~ party + del_gender  + del_age + del_pty_any, data = dat2[dat2$del_gender != "",],
                     family = binomial("logit")))

summary(vet.6 <- glm(veteran ~ party + del_gender  + del_age + del_membership, data = dat2[dat2$del_gender != "",],
                     family = binomial("logit")))

summary(vet.7 <- glm(veteran ~ party + del_gender + del_education + del_age + del_exp_dichotom + del_pty_any + del_membership, data = dat2[dat2$del_education != "" & dat2$del_gender != "",],
                     family = binomial("logit")))


stargazer(list(vet.0, vet.1,  vet.3, vet.4,  vet.6, vet.7),omit = "party",
          covariate.labels = c("Female","Age","University Degree","Legislative Experience",
                               "Party Office (regional/federal)","Length of Party Membership"))

###############################################################################
#                                 Appendix                                    #
###############################################################################


###############################################################################
#                                 Figure 5                                    #
###############################################################################

mm.all <- mm(dat[dat$party_id != "",], all, id = ~ResponseID, alpha = 0.05)

mm.all.w <- mm(dat[dat$party_id != "",], all, id = ~ResponseID, weights = ~freqw, alpha = 0.05)


mm.cdu <- mm(dat[dat$party_id == "CDU",], all, id = ~ResponseID, alpha = 0.05)

mm.spd <- mm(dat[dat$party_id == "SPD",], all, id = ~ResponseID, alpha = 0.05)

mm.fdp <- mm(dat[dat$party_id == "FDP",], all, id = ~ResponseID, alpha = 0.05)

mm.pds <- mm(dat[dat$party_id == "Die LINKE",], all, id = ~ResponseID, alpha = 0.05)

mm.b90 <- mm(dat[dat$party_id == "B�ndnis 90/Die Gr�nen",], all, id = ~ResponseID, alpha = 0.05)

ggData <- data.frame(rbind(cbind(party = "all, n = 296", mm.all, weights = "no Weights"),
                           cbind(party = "all, n = 296", mm.all.w, weights = "with Weights"),
                           cbind(party = "CDU, n = 56", mm.cdu, weights = "no"),
                           cbind(party = "SPD, n = 63", mm.spd, weights = "no"),
                           cbind(party = "FDP, n = 27", mm.fdp, weights = "no"),
                           cbind(party = "LINKE, n = 62", mm.pds, weights = "no"),
                           cbind(party = "B90/Gr�ne, n = 88", mm.b90, weights = "no")))

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

ggData$feature <- as.character(ggData$feature)

ggData$level <- factor(as.character(ggData$level), levels = rev(c("30 Years Old","45 Years Old","65 Years Old",
                                                                  "Female","Male", "Secondary School","High School Diploma",
                                                                  "Technical Profession","Legal Profession","Teacher","Farmer",
                                                                  "District Nomination","No District Nomination",
                                                                  "Member of State Legislature","Municipal Council","None",
                                                                  "7 Months","2 Years", "6 Years",
                                                                  "Leftist Member of Party","Rightist Member of Party",
                                                                  "Centrist Member of Party", "No Office","Party Chair, Village",
                                                                  "Party Chair, Municipality","Member of Regional Board",
                                                                  "Rarely Active","Often Active")))

# plotting
p = ggplot(ggData[ggData$weights != "with Weights",])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_line(data = data.frame(x = 1:29, y = rep(0.5,29)), aes(x = x, y = y), linetype = 2)
p = p + facet_wrap(~ party, nrow =2) + ylab("Marginal Means") + xlab("")
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"))   
p = p + ylab("Marginal Means") + xlab("") 
p


###############################################################################
#                                 Figure 6                                    #
###############################################################################


ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "Marginal Means"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0


p = ggplot(ggData[ggData$party == "all, n = 296",])
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid(.~ weights, scale = "free") + ylab("Marginal Means") + xlab("")
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"))   
p

###############################################################################
#                                 Figure 7                                    #
###############################################################################

mm.first <- mm(dat[!is.na(dat$veteran) & dat$veteran == 0,], all, id = ~ResponseID, alpha = 0.05)

mm.veteran <- mm(dat[!is.na(dat$veteran) & dat$veteran == 1,], all, id = ~ResponseID, alpha = 0.05)

ggData <- data.frame(rbind(cbind(fold = "First-time,\n n = 134", mm.first),
                           cbind(fold = "Veteran,\n n = 116", mm.veteran)))

ggData$level <- factor(as.character(ggData$level), levels = rev(c("30 Years Old","45 Years Old","65 Years Old",
                                                                  "Female","Male", "Secondary School","High School Diploma",
                                                                  "Technical Profession","Legal Profession","Teacher","Farmer",
                                                                  "District Nomination","No District Nomination",
                                                                  "Member of State Legislature","Municipal Council","None",
                                                                  "7 Months","2 Years", "6 Years",
                                                                  "Leftist Member of Party","Rightist Member of Party",
                                                                  "Centrist Member of Party", "No Office","Party Chair, Village",
                                                                  "Party Chair, Municipality","Member of Regional Board",
                                                                  "Rarely Active","Often Active")))
ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

p = ggplot(ggData)
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid( ~ fold, scales = "free") + xlab("") + ylab("")#+ ylab("Marginal Means") 
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                             axis.ticks = element_line(colour = "black"),
                             legend.key = element_rect(fill = "white"),
                             plot.margin = unit(c(0, 0, 0, 0), "cm"),
                             panel.background = element_rect(fill = "white", colour = NA),
                             panel.border = element_rect(fill = NA, colour = "grey50"),
                             panel.grid.major = element_line(colour = "grey90", size = 0.2),
                             panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                             strip.background = element_rect(fill = "white", colour = "grey50"),
                             strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                             strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p


###############################################################################
#                                 Figure 8                                    #
###############################################################################


mm.woman <- mm(dat[dat$del_gender == "Female",], all, id = ~ResponseID)
mm.man <- mm(dat[dat$del_gender == "Male",], all, id = ~ResponseID)
mm.old <- mm(dat[!is.na(dat$del_age) & dat$del_age >= (mean(dat$del_age,na.rm = TRUE)),], all, id = ~ResponseID)
mm.young <- mm(dat[!is.na(dat$del_age) & dat$del_age <= (mean(dat$del_age,na.rm = TRUE)),], all, id = ~ResponseID)

ggData <- data.frame(rbind(cbind(fold = "   Female,    \n n = 119", mm.woman),
                           cbind(fold = "     Male,    \n n = 131", mm.man),
                           cbind(fold = "       Old,      \n n = 104", mm.old),
                           cbind(fold = "    Young,    \n n = 116", mm.young)))

ggData$level <- factor(as.character(ggData$level), levels = rev(c("30 Years Old","45 Years Old","65 Years Old",
                                                                  "Female","Male", "Secondary School","High School Diploma",
                                                                  "Technical Profession","Legal Profession","Teacher","Farmer",
                                                                  "District Nomination","No District Nomination",
                                                                  "Member of State Legislature","Municipal Council","None",
                                                                  "7 Months","2 Years", "6 Years",
                                                                  "Leftist Member of Party","Rightist Member of Party",
                                                                  "Centrist Member of Party", "No Office","Party Chair, Village",
                                                                  "Party Chair, Municipality","Member of Regional Board",
                                                                  "Rarely Active","Often Active")))

ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0


p = ggplot(ggData)
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid( ~ fold, scales = "free") + xlab("") + ylab("")
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                                    axis.ticks = element_line(colour = "black"),
                                    legend.key = element_rect(fill = "white"),
                                    plot.margin = unit(c(0, 0, 0, 0), "cm"),
                                    panel.background = element_rect(fill = "white", colour = NA),
                                    panel.border = element_rect(fill = NA, colour = "grey50"),
                                    panel.grid.major = element_line(colour = "grey90", size = 0.2),
                                    panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                                    strip.background = element_rect(fill = "white", colour = "grey50"),
                                    strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                                    strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p


###############################################################################
#                                 Figure 9                                    #
###############################################################################


mm.high <- mm(dat[dat$del_education == "studied",], all, id = ~ResponseID)
mm.low <- mm(dat[dat$del_education == "not studied",], all, id = ~ResponseID)

mm.right <- mm(dat[!is.na(dat$del_ideology) & dat$del_ideology > 0,], all, id = ~ResponseID)
mm.left <- mm(dat[!is.na(dat$del_ideology) & dat$del_ideology < 0,], all, id = ~ResponseID)

ggData <- data.frame(rbind(cbind(fold = "High-Educated,\n n = 124", mm.high),
                           cbind(fold = "Low-Educated, \n n = 91", mm.low),
                           cbind(fold = "Right-leaning,\n n = 78", mm.right),
                           cbind(fold = "Left-leaning, \n n = 105", mm.left)))

ggData$level <- factor(as.character(ggData$level), levels = rev(c("30 Years Old","45 Years Old","65 Years Old",
                                                                  "Female","Male", "Secondary School","High School Diploma",
                                                                  "Technical Profession","Legal Profession","Teacher","Farmer",
                                                                  "District Nomination","No District Nomination",
                                                                  "Member of State Legislature","Municipal Council","None",
                                                                  "7 Months","2 Years", "6 Years",
                                                                  "Leftist Member of Party","Rightist Member of Party",
                                                                  "Centrist Member of Party", "No Office","Party Chair, Village",
                                                                  "Party Chair, Municipality","Member of Regional Board",
                                                                  "Rarely Active","Often Active")))

ggData$stat_label <- NA
ggData$stat_label[ggData$statistic == "mm"] <- "MM"

ggData$hline <- NA
ggData$hline[ggData$statistic == "mm"] <- 0.5

ggData[is.na(ggData$lower),"lower"] <- 0
ggData[is.na(ggData$upper),"upper"] <- 0

p = ggplot(ggData)
p = p + geom_pointrange(aes(x = as.factor(level), y = estimate, ymin = lower, ymax = upper))
p = p + geom_hline(aes(yintercept = hline), linetype = 2)
p = p + facet_grid( ~ fold, scales = "free") + xlab("") + ylab("")#+ ylab("Marginal Means") 
p = p + coord_flip() + theme(axis.text = element_text(size = rel(0.8)),
                                       axis.ticks = element_line(colour = "black"),
                                       legend.key = element_rect(fill = "white"),
                                       plot.margin = unit(c(0, 0, 0, 0), "cm"),
                                       panel.background = element_rect(fill = "white", colour = NA),
                                       panel.border = element_rect(fill = NA, colour = "grey50"),
                                       panel.grid.major = element_line(colour = "grey90", size = 0.2),
                                       panel.grid.minor = element_line(colour = "grey98", size = 0.5),
                                       strip.background = element_rect(fill = "white", colour = "grey50"),
                                       strip.text.x = element_text(size = 12, colour = "black", face = "bold"),
                                       strip.text.y = element_text(size = 10, colour = "black", face = "bold"))   
p





###############################################################################
#                                Utilities                                    #
###############################################################################

utility <- function(n.sim = 1000, seed = 1104, dat.tmp){
  require(MASS)
  tmp <- lm(all, data = dat.tmp)
  set.seed(seed)
  tmp.coefs <- mvrnorm(n.sim, coef(tmp), vcov(tmp))
  
  age.coefs <- tmp.coefs[,grepl("Old",colnames(tmp.coefs))] # Age
  gender.coefs <- tmp.coefs[,grepl("gender",colnames(tmp.coefs))] # Gender
  education.coefs <- tmp.coefs[,grepl("education",colnames(tmp.coefs))] # Education
  occupation.coefs <- tmp.coefs[,grepl("occupation",colnames(tmp.coefs))] # Occupation
  district.coefs <- tmp.coefs[,grepl("district",colnames(tmp.coefs))] # District nomination
  experience.coefs <- tmp.coefs[,grepl("experience",colnames(tmp.coefs))] # Experience
  length.coefs <- tmp.coefs[,grepl("party_length",colnames(tmp.coefs))] # Party Membership
  office.coefs <- tmp.coefs[,grepl("party_office",colnames(tmp.coefs))] # Party Office
  activity.coefs <- tmp.coefs[,grepl("activity",colnames(tmp.coefs))] # Activity
  ideology.coefs <- tmp.coefs[,grepl("ideology",colnames(tmp.coefs))] 
  
  # choose feature level with largest variance
  age.coefs <- apply(age.coefs, 1, function(x) abs(range(x)[1]- range(x)[2])   )
  occupation.coefs <- apply(occupation.coefs, 1, function(x) abs(range(x)[1]- range(x)[2])   )
  experience.coefs <- apply(experience.coefs, 1, function(x) abs(range(x)[1]- range(x)[2])   )
  office.coefs <- apply(office.coefs, 1, function(x) abs(range(x)[1]- range(x)[2])   )
  length.coefs <- apply(length.coefs, 1, function(x) abs(range(x)[1]- range(x)[2])   )
  ideology.coefs <- apply(ideology.coefs, 1, function(x) abs(range(x)[1]- range(x)[2])   )
  
  sum.coefs <- apply(data.frame(cbind(age.coefs, gender.coefs, education.coefs, occupation.coefs,
                         district.coefs, experience.coefs, length.coefs, office.coefs,
                         ideology.coefs, activity.coefs)), 1, function(x) sum(abs(x)))
  
  importance <- apply(data.frame(cbind(age.coefs, gender.coefs, education.coefs, occupation.coefs,
                         district.coefs, experience.coefs, length.coefs, office.coefs,
                         ideology.coefs, activity.coefs)), 2, function(x) (abs(x)/sum.coefs)*100)

  
  tmp.out <- apply(importance, 2, quantile, probs = c(0.025, 0.5, 0.975))
  
  
  out <- data.frame(rbind(cbind(feature = "Age", importance = tmp.out[2,1], lower = tmp.out[1,1], upper = tmp.out[3,1]),
  cbind(feature = "Gender", importance = tmp.out[2,2], lower = tmp.out[1,2], upper = tmp.out[3,2]),
  cbind(feature = "Education", importance = tmp.out[2,3], lower = tmp.out[1,3], upper = tmp.out[3,3]),
  cbind(feature = "Occupation", importance = tmp.out[2,4], lower = tmp.out[1,4], upper = tmp.out[3,4]),
  cbind(feature = "Nomination", importance = tmp.out[2,5], lower = tmp.out[1,5], upper = tmp.out[3,5]),
  cbind(feature = "Experience", importance = tmp.out[2,6], lower = tmp.out[1,6], upper = tmp.out[3,6]),
  cbind(feature = "Membership", importance = tmp.out[2,7], lower = tmp.out[1,7], upper = tmp.out[3,7]),
  cbind(feature = "Party Office", importance = tmp.out[2,8], lower = tmp.out[1,8], upper = tmp.out[3,8]),
  cbind(feature = "Ideology", importance = tmp.out[2,9], lower = tmp.out[1,9], upper = tmp.out[3,9]),
  cbind(feature = "Activity", importance = tmp.out[2,10], lower = tmp.out[1,10], upper = tmp.out[3,10])))
  
  out[,2] <- as.numeric(as.character(out[,2]))
  out[,3] <- as.numeric(as.character(out[,3]))
  out[,4] <- as.numeric(as.character(out[,4]))
  
  return(out)
}  


u.veteran <- utility(dat.tmp = dat[dat$veteran == 1,])
u.first <- utility(dat.tmp = dat[dat$veteran == 0,])

